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A PRECONDITIONING METHOD FOR SHAPE OPTIMIZATION GOVERNED BY THE 

EULER EQUATIONS * 

EYAL A RIAN * AND VEER N. VATSA * 


Abstract. We consider a classical aerodynamic shape optimization problem subject to the compressible 
Euler flow equations. The gradient of the cost functional with respect to the shape variables is derived with 
the adjoint method at the continuous level. The Hessian (second order derivative of the cost functional with 
respect to the shape variables) is approximated also at the continuous level, as first introduced by Arian 
and Ta’asan (1996). The approximation of the Hessian is used to approximate the Newton step which is 
essential to accelerate the numerical solution of the optimization problem. The design space is discretized 
in the maximum dimension, i.e., the location of each point on the intersection of the computational mesh 
with the airfoil is taken to be an independent design variable. We give numerical examples for 86 design 
variables in two different flow speeds and achieve an order of magnitude reduction in the cost functional at 
a computational effort of a full solution of the analysis partial differential equation (PDE). 

Key words, airfoil, Euler, Hessian, optimal shape, preconditioning. 

Subject classification. Applied and Numerical Mathematics 

1. Introduction. In the last decade there is a growing interest in optimal shape design in conjunction 
with fluid dynamics. The new interest in this classical field is due to advances in computer performance 
and improvements in algorithms for the numerical solution of the flow equations, which makes it practical 
to employ computational fluid dynamics (CFD) codes for shape optimization. 

Historically, the inverse problem of designing a two-dimensional aerodynamical shape profile to attain a 
desired pressure distribution was first studied by Lighthill [2]. However, in general, the desired solution of 
the flow equations on the wing profile is not known in advance and a more difficult optimization problem 
has to be defined. The resulting optimization problems are extremely computationally expensive since it is 
necessary to re-solve the flow equations in each optimization step. Derivative based optimal shape methods 
for problems governed by flow were studied by Pironneau [3] where the adjoint method was used for the 
minimum- drag problem in Stokes flow and subsequently in flows governed by the incompressible Navier- 
Stokes equations. Jameson [4] applied the adjoint method for optimal shape design of an airfoil using the 
Euler equations. Since then there were numerous publications on shape optimization using CFD, most of 
which are contained, or referred to, in Mandel et al. [5] and in the proceedings of the AIAA and ECCOMAS 
meetings. Recently, Soemarwoto [6] computed optimal shapes of airfoils using the compressible Navier- Stokes 
equations, and Jameson et al. [7] obtained results for a three dimensional wing. 

A major effort in that field is in the reduction of the enormous computational effort required to solve 
aerodynamic optimization problems. There are mainly two reasons for the high computational cost; one 
reason is the high computational cost required by each cost functional evaluation, and the second reason 
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19480 while the first author was in residence at the Institute for Computer Applications in Science and Engineering (ICASE), 
M/S 403, NASA Langley Research Center, Hampton, VA 23681-0001. 
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* Aerodynamic and Acoustic Methods Branch, M/S 128, NASA Langley Research Center, Hampton, VA 23681-0001 (email: 
v.n.vatsa@larc.nasa.gov). 
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is the combination of an highly ill-conditioned cost functional Hessian and large design space (hundreds of 
design variables). 

We give here a brief summary of the effort done so far to accelerate the numerical solution of such prob- 
lems. Jameson [4] applied multigrid methods to accelerate the solution of the analysis and adjoint equations, 
while performing the optimization steps only on the finest grid; this shortens the cost functional evaluation 
time. Jameson used a smoothing procedure on the gradients to reduce the ill-conditioning difficulty. Ta’asan 
[8] has proposed the “one-shot” multigrid method where the design space is finite-dimensional and the op- 
timization steps are performed on coarse grids only (Beux and Dervieux [9] proposed a different “one-shot” 
method which is done on a single computational grid but the number of control parameters is progressively 
increased performing a multilevel optimization). Ta’asan et al. [10] and Kuruvila et al. [11] applied these 
ideas to airfoil shape optimization using the full- potential equations. Arian [12] and Arian and Ta’asan 
[13], [14] extended the one-shot multigrid method to infinite dimensional design space in which the optimiza- 
tion steps are performed on all scales including the coarsest and the finest grids. However, the multigrid 
one-shot method requires major change in a given code before it can be applied and in many applications 
this is undesirable; one would like to use a method which is simple to include in a given code. This led 
to the development of new methods which are much simpler to apply and still are very effective. Ta’asan 
[15] proposed a single-grid method to accelerate the numerical convergence of boundary value optimization 
problems, known as the “pseudo- time” method. In the “pseudo- time” method the boundary conditions of 
the state and costate equations are maintained feasible, together with the design equation (defined on the 
boundary only), throughout the computation (see also [16]). Iollo et al. [17] applied the pseudo-time method 
to airfoil shape optimization governed by the Euler equations in two dimensions; the optimum is reported 
to be recovered in the computational effort of roughly 4 full solutions of the analysis problem for 4 and for 
8 design variables. 

One way to lower the computational cost is by approximating the cost functional Hessian and using the 
inverse of the Hessian to act on the gradient. The most popular method to approximate the Hessian, and 
its inverse, very effectively is the Quasi-Newton method (see for example Gill et al. [18]). However, Quasi- 
Newton methods are of low- rank and are difficult to apply efficiently for large ill-conditioned problems. 
The alternative to the Quasi-Newton method is a direct construction of the Hessian matrix. Haftka [19] 
introduced a method to compute the Hessian by solving N sensitivity equations, where N is the number of 
design variables, together with the adjoint equation. Jou et al. [20] use an approximation of such a method 
to construct the Hessian matrix for shape optimization of a three dimensional wing governed by the full- 
potential equation coupled with a boundary layer equation. Arian and Ta’asan [1] introduced a local mode 
analysis of the Hessian, H, for inviscid aerodynamic optimization problems and suggest a preconditioning 
method by transforming the Newton equation in Fourier space, Hs = — to the real space resulting in a 
preconditioning ODE for s ( g here denotes the derivative of the cost functional with respect to the design 
variable). The analysis indicates that inviscid aerodynamic optimization problems are highly ill-conditioned 
and their condition number increases as the design space is refined; it is therefore necessary to approximate 
the Hessian, or rather the effect of its inverse on the gradient (s = —H~ 1 g) in order to obtain effective 
convergence for such problems. 

In this paper, the preconditioning method, suggested by Arian and Ta’asan [1] is re-derived for the 
compressible Euler equations in their conservative form in two space dimensions and is tested numerically 
for an optimal shape design problem in two dimensions. The design space is discretized in the maximum 
dimension, i.e., the location of each point on the intersection of the mesh with the shape is taken to be 
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an independent design variable. Thus, we take the continuous approach into its extreme, both in the 
approximations of the derivatives (gradient and Hessian) and in the discretization of the design space (in 
the continuous approach the derivation of the optimality conditions is done in the continuous, infinite- 
dimensional, level followed by discretization). The optimization algorithm is composed of solving numerically 
the Euler and their adjoint equations followed by a solution of the preconditioning ordinary differential 
equation (ODE) on the profile, and updating the profile’s shape with that solution, resulting in a very 
efficient algorithm. We give numerical examples for 86 design variables in two different flow speeds and 
achieve an order of magnitude reduction in the cost functional at a computational effort of a full solution of 
the analysis PDE. 

2. Problem Definition. Let U denote the vector of state variables: 

(2.1) U = ( p , pu, pv , pw, pE) T . 

In terms of the Jacobian matrices, A ~ ( A, B ), defined in appendix A, the Euler equations are given by 
(conservative form) 

, x v(Am = o inn 

2.2) v ' 

u • n = 0 on T, 

where u = (u,u), T denotes the solid wall (airfoil), and with additional appropriate boundary conditions at 
the far-field. 

The minimization problem in the continuum level is defined by 

(2.3) minF(p) = <f (p(0 ~ P* (Of 

ym Jr 

subject to p satisfying the Euler equations in the domain n. Note that the design variable is composed of 
two functions: y u (x) and y L {x) denoting the upper and lower sections of the airfoil. Fig.(l) depicts the 
geometry of the airfoil. 

The pressure, p, depends on the fourth state variable, by the relation (the parameter 7 is set to 1.4) 

(2.4) p = (7 - 1 )pE - ^-yJpq 2 , 
where q is the flow speed, q = y/u 2 + v 2 . 

3. The Adjoint Equations. We derive the first order necessary conditions for a minimum in the 
continuum level. The derivation is done by taking a variation to the design variable y(£), 

(3.1) y" +1 (0 = y n (0 + ea(0 

and computing its effect on the cost functional with the adjoint method: 

(3.2) F" +1 = F” + ed(0ff(0- 
The derivative p(£) is derived in appendix B: 

(3-3) g(0 = ^M(p - p*) 2 + 2ff (p -P*) ~ (f • n) - £ [(^±^)(tT ■ €)(y ■ n)] =0 on T, 

where R denotes the local radius of curvature, n is a unit vector normal to the surface, y is a unit vector in 
the vertical direction, and A is the vector of costate variables, 

^ 4 x A = (Ai, A 2 , A 3 , A 4 ) t 

A = (A 2 ,A 3 ), 
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Fig. 1. The local coordinate system around the airfoil . 


which satisfies the following “adjoint” equations (with additional far-field boundary conditions): 


(3.5) 


/ T §§ + * T f = 0 infl 
A ■ n + 2 (p — p*) = 0 on I\ 


4. Local Mode Analysis of the Hessian. We derive the symbol of the Hessian for the optimization 
problem, defined by Eqs.(2.2)-(2.3), along the same lines as done by Arian and Ta’asan [1], The derivation 
is divided into subsections with the following content. 

In §4.1 we find the relation between the state variables and the error in the shape. To achieve this we 
use the following representation of the error in the shape, “(0, 


(4.1) 


6(0 = «(«!)«>"* 


and the corresponding errors in the state variables, {Uj}j =1 , 

4 

(4-2) Uj = ]T 

fc= 1 


where v k j denotes the f th component of the null space eigenvector V k , and f3 k are parameters. 

In section §4.2 we find the relation between the symbol of the pressure variable and the symbol of the 
error in the shape variable. 

In §4.3 we use the expression derived in §4.2 to compute the symbol of the Hessian. 


4.1. Null Space of the State Equation Jacobian. We define the matrix L by 


(4.3) 


rr7 A dU n dU 

LU = A—+B — . 
ox oy 


Euler equations for fluids satisfy the relation p = p/(e), where e — E— \q 2 , and as a result AU = 0 (see for 
example Hirsch [21]). Therefore the error equations of (4.3) has the following form: 


(4.4) 


LU 


A dU 8U 

A d^ + B d^’ 
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and the symbol of L is given by 


(4.5) 


L({jJi,u> 2) = iw\A + 


The solution of the linearized Euler error equations is spanned by four eigenvectors that correspond to 
four roots of the equation 

(4.6) Det(L) = ^wfu 2 [2 u\u 2 + (7 - l) 7 (w 2 - 2 E)(wf + w|)] = 0. 

The roots are given by u>\ = {0, }, where we use the relation M = q/c with 


(4.7) 


c 2 = 7(7 - 1) 




For each of the four roots there corresponds an eigenvector, Vfc, for which LVk = 0. 

4.1.1. The root u)\ = 0. The two eigenvectors that correspond to the root u)\ = 0 are given by 


(4.8) 


Vl(uJi = 0) = (uii, ^12^ Vi3, t^i4) T — ( — ^7,0,0, 1) T 

Vliwi = 0) = (V21,V22,V23,V24) T = (f » T 0? 0) T . 


We claim that these vectors do not play a role in the determination of the Hessian’s symbol since they do 
not affect the pressure distribution on the boundary. In order to show this quantitatively we examine the 
relation (2.4) between the pressure, p, and the state variables. The error in p is given by 

(4.9) p= (7-l)(yP~«iPh 1 -u 2 pu2 + pE}. 

For u)\ = 0 the contribution to the error in p by the first two eigenvectors is given by 
2 2 

(4.10) (7 - 1) [(y v n - u i v i2 - U2V13 + 1>!4^ + (7^21 - U1V22 - U2V23 + V24)]a(wi,w 2 )e‘‘‘ ;2T? . 

However, it is easy to check that 


(4.11) 


\vn — U\V \ 2 - U 2 V\3 + V14, = 0 
2 

9 V21 ~ U\V22 ~ ^2^23 H~ ^24 = 0. 


4.1.2. The root uj\ = — • The third eigenvector, which corresponds to the root lj 2 = —u>iV M 2 — 1, 

is not physical since it corresponds to a solution which is diverging in the subsonic regime and propagating 
into the — x direction in the supersonic regime: e xu}2T1 — e _tu;i v/M2 ~ l77 (for subsonic flow, M < 1, the above 
exponent diverges as tj goes to the far field (for more details see Arian and Ta’asan [1])). 

4.1.3. The root We conclude that the only eigenvector that is physical and relevant to 

the Hessian’s symbol is the fourth eigenvector that corresponds to the root u)2 = u\\/M 2 — 1: 

T 

(4.12) V4 = (i>4i,v 4 2, V43,vm) t = ,2(7- 1)(M 2 - 1), -2(7- 1)\/M 2 - 1, deny^j 


where deny = 


2 + (M 2 -2)( 7 -l) 


<?* 
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4.2. The Symbol of the Pressure. By the boundary condition on the solid wall 
(4-13) u-n+a(^ ■ n) - ^(u-£)(n-y) = 0. 

Therefore the following holds for the coefficient 04 (see Eq. (4.2)): 


(4.14) 

By Eq. (4.9) 


= 


{iuiq{y-n) - gf • n)a 

Viz 


(4.15) p = (7 - 1) (y v 41 - qv i2 + v 4 4^/3 4 = q ^ — ~^= — - ct . 

4.3. Derivation of the Hessian’s Symbol. The cost functional (2.3) can be presented as an inner 
product in L 2 (T): 


(4.16) F = (p-p m ,p-p*) L 2 { r) . 

In order to estimate the Hessian we need to expand the cost functional, F, as a Taylor series in a perturbation 
of the shape, d. A perturbation of the shape, a, results in a perturbed cost functional of the form 

(4.17) F(c* + q:) = /((p-p*)+p+a|?),((p — p *)+p + Q^)(i + ^)\ + h.o.t. 

\ ay dy R f l 2 {V) 

where R is local radius of curvature. 

Let us assume that the operator Z maps the change in the shape, d, with the corresponding change in 
the pressure, p : 


(4.18) p = Za. 

In terms of Z the Taylor expansion of (4.16) is given by: 

F{a + a) = F(a) + ((p - p m )(Z + Z* + 2§£), a) 


lh n 


+ ( [ (z * + % ){z + + h {p ~ p * ){2 ^ + + z) h & ) L H d + hoL 


(4.19) ' \ L v- ' dy >y- ■ Q y / ■ R >y- dy 

Therefore, in terms of Z the symbol of the Hessian is approximated by 

dP,** . ^ , , d P^2 , 1 /_ _*', n dp 


z-z + f y (Z- + Z) + Cf/ + ^P -Pl(2f y + z- + i)). 


% ■n-iu 1 q(y-n) 


(4.20) H = 2( 

where by Eq.(4.15) the symbol of Z is equal to 

^ 4 ' 21 ^ ~ ’ s/W - 1 

We conclude that the symbol of the Hessian is approximated by 

(4.22) H = aJ 1 + 6, 

where a and b are given by 

„ _ o/ V(g *) a \ 

b = 2 [(Re(Z) + g) 3 + £(p - p*)(g + *e(i))). 

TVansforming the Newton equation, Hs = — g, back into the real space we obtain the following precondi- 
tioning ODE, which approximates the Newton step, on the designed profile: 


(4.23) 


(4.24) 


1 S 


+ bs == — g on T. 
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5. The Optimization Algorithm. Starting with an approximation ((7 n (r n ), A n (r n )) (from here on 
the superscript n will be omitted) 

1. Solve partially for U n + l the Euler equations (2.2) (reduce residuals by two 
orders of magnitude) . 

2. Solve partially for A n+1 the adjoint equations (3 . 5) , (reduce residuals by two 
orders of magnitude) . 

3. Compute the gradient, g , with Eq.(3.3). 

4. Solve the Preconditioning equation (4.24) for s. 

5. Update the shape with P 16 ™ = T old + Ss and re-mesh, where 6 is a parameter. 

6. If the L 2 norm of the gradient, ||<7||i, 2 (r), is smaller than some prescribed value 
then stop, else goto 1. 


6. Numerical Test. A detailed description of the numerical scheme that we used for the solution 
of the Euler equations is given in Vatsa and Wedan [22] and in Vatsa et al. [23]. Here we give only a 
brief summary (the adjoint system (3.5) was solved in a similar manner to the solution of Euler system). 
We used a semi-discrete finite-volume algorithm based on Jameson’s Runge-Kutta time-stepping scheme 
(Jameson et al. [24] and Jameson and Baker [25]) to obtain the steady-state of the governing Euler and their 
adjoint equations. A controlled amount of artificial dissipation was added to the central-difference scheme 
to suppress odd-even point decoupling and oscillations in the vicinity of shock waves and stagnation points. 
The discrete design space was composed of the coordinates of the grid intersection with the profile, {yi,i}£Li, 
with N = 89 (Fig.(l) depicts the grid around the airfoil). The leading edge and two trailing edge grid points 
were fixed during the computation and therefore the actual number of design variables was 86. We minimize 
the pressure matching cost functional, (2.3), where the discrete target pressure distribution, {p*}^, was 
computed on the NACA 0012 airfoil, and the initial geometry was taken to be the Onera M6 airfoil. 

Our assumption is that in the advanced phase of aircraft design the wing’s planform is fixed and therefore 
we fixed the leading and trailing edges of the airfoil. The airfoil was divided into two sections: the upper 
part y u (x) and the lower part y L (x). On each section we used 43 design variables, having a total number of 
86 design variables, and the optimization algorithm presented in Sec. 6 was applied. 

Given the discretized gradient along the airfoil, we solve the discrete version of Eq.(4.24) for 

{ 5 i l }il 1 on the airfoil, with the boundary conditions s h = 0 on the trailing and leading edges (our numerical 
experiments show that setting b = 0 in Eq.(4,24) results in a more robust preconditioner, although less 
accurate in some cases). Thus Eq.(4.24) is solved twice, i.e., once on the upper section of the airfoil and 
once on the lower section. At transonic flow, where there are sonic points along the airfoil we introduce a 
cutoff on the denominator of Eq.(4.24): |1 — M 2 \ > e (we used e = 0.4). The solution of Eq.(4.24), {sf }£L lT 
is the new “preconditioned” search direction. 

Two flow speeds were tested: subsonic with M = 0.5 and transonic with M = 0.75. 

6.1. Subsonic Flow. Table 6.1 depicts the cost functional history in logarithmic scale for Mach number 
0.5. At each optimization step we used a fixed number of iterations (50) for the flow and adjoint partial 
solutions. We observe an average reduction of about 0.23 in the cost functional at every iteration, while 
the computational effort required by each iteration is of about a third of one cost functional evaluation (50 
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iterations for the state equations and another 50 iterations for the costate equations, while it takes about 
300 iterations to solve the Euler equations to machine zero accuracy (one cost functional evaluation)). 

Fig.2 depicts the error (err), preconditioned solution (s), and gradient (g), for the first four iterations 
(all normalized). The preconditioning equation (4.24) transformed the gradient (g) into a new function 
(s) which approximates the actual error in the shape (err). The actual error in the shape is computed by 
{err± = y* ar 9 €t _ where y t . ar9et denotes the y coordinates at the mesh point i along the boundary 

of the NACA 0012 airfoil, and denotes the y coordinates of the current iterate airfoil (starting with the 
Onera M6 airfoil) . 


Iteration 

1 

2 

3 

4 

5 

Cost Functional 

4.32 10" 4 

1.17 10“ 4 

2.16 1(T 5 

6.04 10" 6 

1.28 HT 6 


Table 6.1 ” ” ~~ 

Cost functional history for Mach number 0 . 5 . The initial airfoil is Onera M6 and the target pressure distribution is taken 
from a NACA 0012 airfoil. 


6.2. Transonic Flow. The results for the transonic case are depicted in table 6.2 and Fig.3. 

Table 6.2 depicts the cost functional history in logarithmic scale for Mach number 0.75. (The initial 
airfoil is Onera M6 and the target pressure distribution is taken from a 0012 airfoil.) At each optimization 
step we used a fixed number of iterations (100) for the flow and adjoint partial solutions. In this case the 
convergence rate is slower than in the subsonic flow; we observe an average reduction of about 0.33 in the 
cost functional at every iteration, while the computational effort required by each iteration is of about a two 
thirds of one cost functional evaluation. 

Fig.3 depicts the error (err), preconditioned (s), and gradient (g), for the first four iterations (all nor- 
malized). 


Iteration 

. 1 

2 

3 

4 

5 

Cost Functional 

6.52 10~ 3 

2.04 10~ 3 

5.04 10 -4 

2.20 10 -4 

7.64 10“ 5 


Table 6.2 

Cost functional history for Mach number 0 . 75 . The initial airfoil is Onera M6 and the target pressure distribution is 
taken from a NACA 0012 airfoil. 


7. Discussion and Concluding Remarks. An important issue in the field of numerical methods for 
optimization problems governed by PDEs is the discrete versus the continuum approaches. In the discrete 
approach, the derivation of the first and second order sensitivities is done in the discrete level (for example 
with automatic differentiation), while in the continuum approach the derivation is done in the continuous, 
infinite-dimensional, level followed by discretization. In this paper we take the continuous approach to its 
extreme; both in the approximations of the derivatives (gradient and Hessian) and in the discretization of the 
design space. The main advantage of the continuous approach is that we can use the structure of the problem 
by analyzing it in the PDE level, resulting in a very efficient computational method, as demonstrated in 
this paper. However, we find that the continuous approach has its disadvantages. We had some difficulties 
with the point-wise control of the shape, and with the sensitivity gradients accuracy, around singular points. 
In our opinion, there is a tradeoff between robustness and efficiency of the numerical scheme, where the 
continuous approach leads to schemes which are less robust but much more computationally efficient than 
the schemes resulting from the discrete approach. 
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Appendix A. Definition of the Euler Jacobian Matrices. 

0 1 0 0 

— u 2 + J Y^q 2 (3 — 7)u — (7— l)v 7— 1 

— uv v u 0 

-u(yE — (7 — 1)<Z 2 ) jE — J ^ L {q 2 + 2u 2 ) — (7— 1 )uv 7 u 

0 0 10 
— uv v u 0 

-v 2j r :L Y L q 2 — (7— l)u (3 — 7)1; 7-I 

~v(^yE - (7 - l)<? 2 j — (7— 1 )uv *yE — J ^ L {q 2 + 2v 2 ) *yv 


Appendix B. Derivation of the Shape Derivatives. We compute the derivative of the cost func- 
tional with respect to the shape of the profile, t/(£), with the adjoint method in the continuous level. A 
Lagrangian is defined as follows: 

(B.l) C(U, A) = j> (p- p*) 2 d£ + J [A • (A^ + B—f^dxdy. 

Then an arbitrary perturbation of the profile’s shape, d, is introduced with a small amplitude e : 

(B.2) y(x) new = y(x) old + ea(x) 

The resulting first order variation of the Lagrangian, is given by (see for example [3], [12] for further details) 

1 = [ 2<# + <5 |> <p “ ■ *)]« + jf [0 ■ + B T ^)]dxdy 

(B.3) + J A T [(n x A + n 2 B)^]de, 
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where n\ and n 2 are the components of the unit normal vector in the x and y directions, respectively, and 
where we have used the following relation (see [3]) 

(B.4) dC ew = + 

Next we compute the term A 7 (nj A 4- n 2 B^jU on the solid wall. The no-slip boundary condition u ■ n = 0 
implies 

(B.5) u- n + d(|f -n) - §|(y-n)(u •£) = 0 on T, 

where we have used the following relation 

—* 

(B.6) = n old -—{y-n)t. 

Using relation (B.5) in U on the solid wall, and the variational form of Eq.(2.4), results in the following 
equation: 

(B.7) A T (n 1 A + n 2 B)^ = (A-n)p-^^((f n)d-f (y-n)(«-0) on T. 

Equating the terms in £ that multiply u on T to zero results in the solid wall boundary condition for 
the adjoint problem: 

(B.8) A • n + 2(p — p*) = 0 on T. 

Equating the terms in £ that multiply a on T to zero results in the design equation: 

(B.9 )g(0 =^(p- P*) 2 + 2 g(p - p*) - (*^) (If . n) - ^ [(M±^)( 5 . 0(y-n)] =0 on T. 
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Fig. 2. The vertical axis represents the normalized actual error (err) in the airfoiTs shape, the normalized estimated error 
(s), with the new preconditioning method, and the normalized sensitivity gradient (g) of a cost functional with respect to the 
shape parameters along the airfoil. The horizontal axis represents the location of the computational mesh points on the airfoil. 
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